То, что мы в быту привыкли называть распределением - это функция плотности вероятности.
Плотность вероятности - это функция, описывающая вероятность получения разных значений случайной величины
\[f(y;\mu, \sigma)= \frac {1}{\sigma \sqrt{2 \pi}} e^{-\frac{(y-\mu)^2}{2\sigma^2}}\]
Среднее:   \(E(Y) = \mu\)
Дисперсия: \(var(Y) = \sigma^2\)
\(-\infty \le Y \le +\infty\)
\[f(y;\mu)= \frac{\mu^y \times e{-\mu}}{y!}\]
Среднее:   \(E(Y) = \mu\)
Дисперсия: \(var(Y) = \mu\)
Важное свойство: При увеличении значения \(\mu\) увеличивается размах варьирования
\(0 \le Y \le +\infty\),
\(Y\) целочисленные
\(f(y; \mu, \nu) = \frac{1}{\Gamma(\nu)}\times (\frac{\nu}{\mu})^{\nu} \times y^{\nu-1} \times e^{\frac{y \times \nu}{\mu}}\)
Среднее:   \(E(Y) = \mu\)
Дисперсия: \(var(Y) = \frac {\mu^2}{\nu}\)
Параметр \(\nu\) определяет степень избыточности дисперсии
\(0 < Y \le +\infty\)
Внимание! \(Y\) строго больше 0
\(f(y; k, \mu) = \frac{\Gamma(y + k)}{\Gamma(k) \times \Gamma(y+1)} \times (\frac{k}{\mu + k})^k \times (1 - \frac{k}{\mu + k})^y\)
Это смесь Пуассоновского и Гамма распределений: \(Y\) демонстрируют распределение Пуассона с \(\mu\), подчиняющимися Гамма-распределению.
Среднее:   \(E(Y) = \mu\)
Дисперсия: \(var(Y) = \mu + \frac {\mu^2}{k}\)
Параметр \(k\) определяет степень избыточности дисперсии.
Важное свойство: Приближается к распр. Пуассона при очень больших \(k\).
\(0 \le Y \le +\infty\),   \(Y\) целочисленные
\(f(y; N, \pi) = \frac{N!}{y! \times (N-y)!} \times \pi^y \times (1 - \pi)^{N-y}\)
Среднее:    \(E(Y) = N \times \pi\)
Дисперсия: \(var(Y) = N \times \pi \times (1-\pi)\)
Параметр \(N\) определяет количество объектов в испытании
Парметр \(\pi\) - вероятность события (\(y = 1\))
\(0 \le Y \le +\infty\),   \(Y\) целочисленные
Предсказаные значения, т.е. \(E(Y_i)\), лежат на линии регрессии
из кн. Zuur et al.2009
Для каждого типа распределения характер связи между \(\mu_i\) и значениями предиктора может быть разным.
Функция, описывающая связь между \(\mu_i\) и значениями предикторов, называется связующей функцией
| Характер величин | Распределение | Связующая функция (link function) |
|---|---|---|
| Непрерывные величины, потенциально варьирующие в пределах \(-\infty , + \infty\) | Гауссовское (Нормальное) | identity \(X\beta = \mu\) |
| Бинарные величины (1; 0), или количество (доля) объектов, относящихся к одному из двух классов | Биномиальное распределение | logit \(X\beta = ln(\frac{\mu}{1 - \mu})\) |
| Счетные величины (0, 1, 2, 3…) | Распределение Пуассона или Отрицательное биномиальное распределение | log \(X\beta = ln(\mu)\) |
M_sim1 <- glm(y ~ x, data = dat, family = poisson(link = "log"))
Даные взяты из работы Khaitov, 2013
juv_ad <- read.table("data/mussel_juv_ad.csv", sep=";", header =T)
head(juv_ad, 12)
## Year Bank Sample Juv Adult ## 1 1997 vor4 1 90 42 ## 2 1997 vor4 2 134 25 ## 3 1997 vor4 3 166 27 ## 4 1997 vor4 4 168 48 ## 5 1997 vor4 5 102 25 ## 6 1997 vor4 6 69 57 ## 7 1998 vor4 1 410 21 ## 8 1998 vor4 2 73 76 ## 9 1998 vor4 3 347 0 ## 10 1998 vor4 5 402 5 ## 11 1998 vor4 6 158 28 ## 12 1999 vor4 1 282 0
Как можно учесть в модели многолетний характер сбора материала?
M1 <- glm(Juv ~ Adult * factor(Year), data = juv_ad) drop1(M1, test = "F")
## Single term deletions ## ## Model: ## Juv ~ Adult * factor(Year) ## Df Deviance AIC F value Pr(>F) ## <none> 256799 1024 ## Adult:factor(Year) 14 361458 1026 1.72 0.076 . ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M2 <- lm(Juv ~ Adult + factor(Year), data = juv_ad) library(car) Anova(M2)
## Anova Table (Type II tests) ## ## Response: Juv ## Sum Sq Df F value Pr(>F) ## Adult 63886 1 12.9 0.00059 *** ## factor(Year) 1357772 14 19.6 < 2e-16 *** ## Residuals 361458 73 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M2_diag <- fortify(M2) ggplot(M2_diag, aes(x = .fitted, y = .stdresid)) + geom_point() + geom_hline(yintercept = 0) + geom_smooth(se = FALSE)
M3 <- glm(Juv ~ Adult * factor(Year),
data = juv_ad, family = "poisson")
library(car)
Anova(M3)
## Analysis of Deviance Table (Type II tests) ## ## Response: Juv ## LR Chisq Df Pr(>Chisq) ## Adult 272 1 <2e-16 *** ## factor(Year) 5031 14 <2e-16 *** ## Adult:factor(Year) 439 14 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M3_diag <- data.frame(.fitted = predict(M3),
.pears_resid = residuals(M3, type = "pearson"))
Рассеяние остатков выглядит лучше!
В Пуассоновской регрессии мы моделируем изменение распределения Пуассона в зависимости от каких-то предикторов.
В распределении Пуассона \(E(Y) = \mu\) и \(var(Y) = \mu\)
Если в имеющихся данных \(var(Y) > \mu\), то нарушается условие применимости пуассоновской регрессии.
Resid_M3 <- resid(M3, type = "pearson") # Пирсоновские остатки N <- nrow(juv_ad) # Объем выборки p <- length(coef(M3)) # Число параметров в модели df <- (N - p) # число степенейсвободы fi <- sum(Resid_M3^2) /df #Величина fi показывает во сколько раз в среднем sigma > mu для данной модели fi
## [1] 16.8
Дисперсия в 16.756 раза больше среднего!
library(qcc) qcc.overdispersion.test(juv_ad$Juv, type = "poisson")
## ## Overdispersion test Obs.Var/Theor.Var Statistic p-value ## poisson data 103 9079 0
summary(M3)
## ## Call: ## glm(formula = Juv ~ Adult * factor(Year), family = "poisson", ## data = juv_ad) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -10.779 -2.370 -0.074 2.127 7.858 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 5.14952 0.11551 44.58 < 2e-16 *** ## Adult -0.00955 0.00305 -3.13 0.00173 ** ## factor(Year)1998 0.89375 0.12042 7.42 1.2e-13 *** ## factor(Year)1999 0.46774 0.12279 3.81 0.00014 *** ## factor(Year)2000 -0.67585 0.20462 -3.30 0.00096 *** ## factor(Year)2001 0.09944 0.23147 0.43 0.66748 ## factor(Year)2002 -0.46104 0.25068 -1.84 0.06589 . ## factor(Year)2003 0.99365 0.12334 8.06 7.9e-16 *** ## factor(Year)2004 0.34742 0.13293 2.61 0.00896 ** ## factor(Year)2005 0.51396 0.13149 3.91 9.3e-05 *** ## factor(Year)2006 0.40032 0.12997 3.08 0.00207 ** ## factor(Year)2007 -1.51444 0.21138 -7.16 7.8e-13 *** ## factor(Year)2008 0.89713 0.12124 7.40 1.4e-13 *** ## factor(Year)2009 0.47749 0.12401 3.85 0.00012 *** ## factor(Year)2010 0.31445 0.13054 2.41 0.01600 * ## factor(Year)2011 1.05365 0.12258 8.60 < 2e-16 *** ## Adult:factor(Year)1998 -0.01150 0.00336 -3.42 0.00063 *** ## Adult:factor(Year)1999 0.01496 0.01048 1.43 0.15355 ## Adult:factor(Year)2000 0.02677 0.00580 4.62 3.9e-06 *** ## Adult:factor(Year)2001 -0.00481 0.00490 -0.98 0.32698 ## Adult:factor(Year)2002 0.00534 0.00533 1.00 0.31643 ## Adult:factor(Year)2003 -0.02341 0.00411 -5.69 1.2e-08 *** ## Adult:factor(Year)2004 -0.01631 0.00591 -2.76 0.00575 ** ## Adult:factor(Year)2005 -0.00882 0.00381 -2.31 0.02062 * ## Adult:factor(Year)2006 -0.04550 0.00760 -5.99 2.1e-09 *** ## Adult:factor(Year)2007 0.05039 0.00944 5.34 9.4e-08 *** ## Adult:factor(Year)2008 0.00569 0.00349 1.63 0.10312 ## Adult:factor(Year)2009 -0.00642 0.00508 -1.26 0.20637 ## Adult:factor(Year)2010 0.00445 0.00398 1.12 0.26394 ## Adult:factor(Year)2011 0.01949 0.00351 5.55 2.8e-08 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 7992.3 on 88 degrees of freedom ## Residual deviance: 1008.2 on 59 degrees of freedom ## AIC: 1695 ## ## Number of Fisher Scoring iterations: 4
Если избыточнсть дисперсии мнимая, то ее надо устранить, введя в модель соответствующие поправки.
Если избыточность дисперсии истинная, то необходима более серьезная коррекция модели.
Отличие от пуассоновсой модели заключается лишь в том, что в квази-пуассоновских моделях вводится поправка для связи дисперсии и матожидания.
В этой модели матожидание \(E(Y) = \mu\) и дисперсия \(var(Y) = \phi \times \mu\)
Величина \(\phi\) показывает во сколько раз дисперсия превышает матожидание.
\[\phi = \frac{var(Y)}{\mu}=\frac {\frac{\sum{(\epsilon_i)^2}}{N - p}} {\mu} = \frac{\sum{(\epsilon_{pearson})^2}}{N - p}\]
M4 <- glm(Juv ~ Adult * factor(Year), data = juv_ad, family = "quasipoisson") Anova(M4)
## Analysis of Deviance Table (Type II tests) ## ## Response: Juv ## LR Chisq Df Pr(>Chisq) ## Adult 16.2 1 0.000056 *** ## factor(Year) 300.3 14 < 2e-16 *** ## Adult:factor(Year) 26.2 14 0.024 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Важно: Распределение разности девианс лишь приблизительно описывается распределением \(\chi^2\). Поэтому уровни значимости близкие к 0.05 нельзя рассматривать как надежные.
Уровень значимости для взаимодействия в данной модели близок к 0.05. Можно подумать об упрощении модели!
Упрощенная модель
M4a <- glm(Juv ~ Adult + factor(Year), data = juv_ad, family = "quasipoisson") Anova(M4a)
## Analysis of Deviance Table (Type II tests) ## ## Response: Juv ## LR Chisq Df Pr(>Chisq) ## Adult 13.8 1 0.0002 *** ## factor(Year) 254.9 14 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(MASS) M5 <- glm.nb(Juv ~ Adult*factor(Year) , data = juv_ad, link = "log") Anova(M5)
## Analysis of Deviance Table (Type II tests) ## ## Response: Juv ## LR Chisq Df Pr(>Chisq) ## Adult 27 1 0.00000016 *** ## factor(Year) 370 14 < 2e-16 *** ## Adult:factor(Year) 34 14 0.0023 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Уровень значимости для взаимодействия заметно ниже 0.05. Взаимодействие факторов отбросить нельзя!
Проверьте на избыточность дисперсии модель, основанную на отрицательном биномиальном распеделении
Resid_M5 <- resid(M5, type = "pearson") # Пирсоновские остатки N <- nrow(juv_ad) # Объем выборки p <- length(coef(M5)) +1 # Число параметров в модели df <- (N - p) # число степенейсвободы fi <- sum(Resid_M5^2) /df #Величина fi показывает # во сколько раз в среднем sigma > mu для данной модели fi
## [1] 1.47
M5_diag <- data.frame(.fitted = predict(M5),
.pears_resid = residuals(M5, type = "pearson"))
ggplot(M5_diag, aes(x=.fitted, y = .pears_resid)) + geom_point() +
geom_hline(yintercept = 0) + geom_smooth(se = F)
Визуализируйте предсказания модели, основанной на отрицательном биномиальном распределении
MyData <- expand.grid(Year = factor(seq(min(juv_ad$Year),
max(juv_ad$Year))), Adult = seq(min(juv_ad$Adult),
max(juv_ad$Adult)))
MyData$Predicted <- predict(M5, newdata = MyData, type = "response")
ggplot(MyData, aes(x = Adult, y = Predicted, group = Year)) +
geom_line(color = "blue") + geom_hline(yintercept = 0) +
ylab("Ожидаемое количество молоди") +
geom_point(data = juv_ad, aes(x = Adult, y = Juv,
color = factor(Year))) + facet_wrap(~Year,
ncol = 5) + guides(color = FALSE)
data(Affairs, package = "AER") af <- Affairs
Оригинальная работа: Fair, R.C. (1978). A Theory of Extramarital Affairs. Journal of Political Economy, 86, 45–61.
affairs - Количество внебрачных свзяей за последний годgender - полage - возрастyearsmarried - сколько ле в бракеchildren - наличие детейreligiousness - уровеь религиозностиeducation - уровень образованияrating - самооценка ощущений от брака
Mod <- glm(affairs ~ gender * yearsmarried * children * religiousness +
education, data = af, family = "poisson")
summary(Mod)
## ## Call: ## glm(formula = affairs ~ gender * yearsmarried * children * religiousness + ## education, family = "poisson", data = af) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -3.27 -1.78 -1.26 -0.60 6.74 ## ## Coefficients: ## Estimate Std. Error ## (Intercept) 0.6798 0.5380 ## gendermale -1.2262 0.7071 ## yearsmarried 0.2747 0.0666 ## childrenyes 0.7479 0.6325 ## religiousness -0.5268 0.1915 ## education -0.0274 0.0149 ## gendermale:yearsmarried 0.0325 0.0796 ## gendermale:childrenyes 1.7851 0.8714 ## yearsmarried:childrenyes -0.2238 0.0740 ## gendermale:religiousness 0.5167 0.2605 ## yearsmarried:religiousness -0.0183 0.0241 ## childrenyes:religiousness 0.1421 0.2405 ## gendermale:yearsmarried:childrenyes -0.0833 0.0902 ## gendermale:yearsmarried:religiousness -0.0282 0.0288 ## gendermale:childrenyes:religiousness -0.7454 0.3216 ## yearsmarried:childrenyes:religiousness 0.0216 0.0266 ## gendermale:yearsmarried:childrenyes:religiousness 0.0499 0.0324 ## z value Pr(>|z|) ## (Intercept) 1.26 0.2064 ## gendermale -1.73 0.0829 ## yearsmarried 4.13 0.000037 ## childrenyes 1.18 0.2371 ## religiousness -2.75 0.0059 ## education -1.84 0.0662 ## gendermale:yearsmarried 0.41 0.6829 ## gendermale:childrenyes 2.05 0.0405 ## yearsmarried:childrenyes -3.03 0.0025 ## gendermale:religiousness 1.98 0.0473 ## yearsmarried:religiousness -0.76 0.4473 ## childrenyes:religiousness 0.59 0.5544 ## gendermale:yearsmarried:childrenyes -0.92 0.3555 ## gendermale:yearsmarried:religiousness -0.98 0.3274 ## gendermale:childrenyes:religiousness -2.32 0.0205 ## yearsmarried:childrenyes:religiousness 0.81 0.4174 ## gendermale:yearsmarried:childrenyes:religiousness 1.54 0.1242 ## ## (Intercept) ## gendermale . ## yearsmarried *** ## childrenyes ## religiousness ** ## education . ## gendermale:yearsmarried ## gendermale:childrenyes * ## yearsmarried:childrenyes ** ## gendermale:religiousness * ## yearsmarried:religiousness ## childrenyes:religiousness ## gendermale:yearsmarried:childrenyes ## gendermale:yearsmarried:religiousness ## gendermale:childrenyes:religiousness * ## yearsmarried:childrenyes:religiousness ## gendermale:yearsmarried:childrenyes:religiousness ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 2925.5 on 600 degrees of freedom ## Residual deviance: 2518.8 on 584 degrees of freedom ## AIC: 3047 ## ## Number of Fisher Scoring iterations: 7
# Проверка на Overdispersion Resid_Mod <- resid(Mod, type = "pearson") N <- nrow(af) p <- length(coef(Mod)) df <- (N - p) fi <- sum(Resid_Mod^2) /df fi
## [1] 6.44
Mod1 <- glm(affairs ~ gender * yearsmarried * children * religiousness *
education, data = af, family = "quasipoisson")
summary(Mod1)
## ## Call: ## glm(formula = affairs ~ gender * yearsmarried * children * religiousness * ## education, family = "quasipoisson", data = af) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -3.356 -1.699 -1.197 -0.518 6.833 ## ## Coefficients: ## Estimate ## (Intercept) -0.72294 ## gendermale 6.32757 ## yearsmarried -0.13878 ## childrenyes 21.26596 ## religiousness 2.68502 ## education 0.10930 ## gendermale:yearsmarried 1.87507 ## gendermale:childrenyes -27.44754 ## yearsmarried:childrenyes -1.24552 ## gendermale:religiousness -3.35695 ## yearsmarried:religiousness -0.01597 ## childrenyes:religiousness -10.10729 ## gendermale:education -0.52222 ## yearsmarried:education 0.01782 ## childrenyes:education -1.44076 ## religiousness:education -0.23255 ## gendermale:yearsmarried:childrenyes -0.86718 ## gendermale:yearsmarried:religiousness -0.47529 ## gendermale:childrenyes:religiousness 11.26910 ## yearsmarried:childrenyes:religiousness 0.56754 ## gendermale:yearsmarried:education -0.09418 ## gendermale:childrenyes:education 1.97936 ## yearsmarried:childrenyes:education 0.08050 ## gendermale:religiousness:education 0.27615 ## yearsmarried:religiousness:education 0.00383 ## childrenyes:religiousness:education 0.71028 ## gendermale:yearsmarried:childrenyes:religiousness -0.00469 ## gendermale:yearsmarried:childrenyes:education 0.01680 ## gendermale:yearsmarried:religiousness:education 0.01994 ## gendermale:childrenyes:religiousness:education -0.82013 ## yearsmarried:childrenyes:religiousness:education -0.04129 ## gendermale:yearsmarried:childrenyes:religiousness:education 0.01511 ## Std. Error ## (Intercept) 11.55629 ## gendermale 15.47884 ## yearsmarried 1.28616 ## childrenyes 14.96782 ## religiousness 4.61169 ## education 0.77941 ## gendermale:yearsmarried 2.14783 ## gendermale:childrenyes 19.12367 ## yearsmarried:childrenyes 1.45496 ## gendermale:religiousness 5.84741 ## yearsmarried:religiousness 0.47583 ## childrenyes:religiousness 5.74652 ## gendermale:education 1.00639 ## yearsmarried:education 0.08122 ## childrenyes:education 1.01421 ## religiousness:education 0.31845 ## gendermale:yearsmarried:childrenyes 2.31968 ## gendermale:yearsmarried:religiousness 0.70212 ## gendermale:childrenyes:religiousness 7.13962 ## yearsmarried:childrenyes:religiousness 0.53301 ## gendermale:yearsmarried:education 0.12702 ## gendermale:childrenyes:education 1.24743 ## yearsmarried:childrenyes:education 0.09344 ## gendermale:religiousness:education 0.38721 ## yearsmarried:religiousness:education 0.03019 ## childrenyes:religiousness:education 0.39169 ## gendermale:yearsmarried:childrenyes:religiousness 0.76524 ## gendermale:yearsmarried:childrenyes:education 0.13881 ## gendermale:yearsmarried:religiousness:education 0.04196 ## gendermale:childrenyes:religiousness:education 0.46874 ## yearsmarried:childrenyes:religiousness:education 0.03419 ## gendermale:yearsmarried:childrenyes:religiousness:education 0.04622 ## t value ## (Intercept) -0.06 ## gendermale 0.41 ## yearsmarried -0.11 ## childrenyes 1.42 ## religiousness 0.58 ## education 0.14 ## gendermale:yearsmarried 0.87 ## gendermale:childrenyes -1.44 ## yearsmarried:childrenyes -0.86 ## gendermale:religiousness -0.57 ## yearsmarried:religiousness -0.03 ## childrenyes:religiousness -1.76 ## gendermale:education -0.52 ## yearsmarried:education 0.22 ## childrenyes:education -1.42 ## religiousness:education -0.73 ## gendermale:yearsmarried:childrenyes -0.37 ## gendermale:yearsmarried:religiousness -0.68 ## gendermale:childrenyes:religiousness 1.58 ## yearsmarried:childrenyes:religiousness 1.06 ## gendermale:yearsmarried:education -0.74 ## gendermale:childrenyes:education 1.59 ## yearsmarried:childrenyes:education 0.86 ## gendermale:religiousness:education 0.71 ## yearsmarried:religiousness:education 0.13 ## childrenyes:religiousness:education 1.81 ## gendermale:yearsmarried:childrenyes:religiousness -0.01 ## gendermale:yearsmarried:childrenyes:education 0.12 ## gendermale:yearsmarried:religiousness:education 0.48 ## gendermale:childrenyes:religiousness:education -1.75 ## yearsmarried:childrenyes:religiousness:education -1.21 ## gendermale:yearsmarried:childrenyes:religiousness:education 0.33 ## Pr(>|t|) ## (Intercept) 0.950 ## gendermale 0.683 ## yearsmarried 0.914 ## childrenyes 0.156 ## religiousness 0.561 ## education 0.889 ## gendermale:yearsmarried 0.383 ## gendermale:childrenyes 0.152 ## yearsmarried:childrenyes 0.392 ## gendermale:religiousness 0.566 ## yearsmarried:religiousness 0.973 ## childrenyes:religiousness 0.079 ## gendermale:education 0.604 ## yearsmarried:education 0.826 ## childrenyes:education 0.156 ## religiousness:education 0.466 ## gendermale:yearsmarried:childrenyes 0.709 ## gendermale:yearsmarried:religiousness 0.499 ## gendermale:childrenyes:religiousness 0.115 ## yearsmarried:childrenyes:religiousness 0.287 ## gendermale:yearsmarried:education 0.459 ## gendermale:childrenyes:education 0.113 ## yearsmarried:childrenyes:education 0.389 ## gendermale:religiousness:education 0.476 ## yearsmarried:religiousness:education 0.899 ## childrenyes:religiousness:education 0.070 ## gendermale:yearsmarried:childrenyes:religiousness 0.995 ## gendermale:yearsmarried:childrenyes:education 0.904 ## gendermale:yearsmarried:religiousness:education 0.635 ## gendermale:childrenyes:religiousness:education 0.081 ## yearsmarried:childrenyes:religiousness:education 0.228 ## gendermale:yearsmarried:childrenyes:religiousness:education 0.744 ## ## (Intercept) ## gendermale ## yearsmarried ## childrenyes ## religiousness ## education ## gendermale:yearsmarried ## gendermale:childrenyes ## yearsmarried:childrenyes ## gendermale:religiousness ## yearsmarried:religiousness ## childrenyes:religiousness . ## gendermale:education ## yearsmarried:education ## childrenyes:education ## religiousness:education ## gendermale:yearsmarried:childrenyes ## gendermale:yearsmarried:religiousness ## gendermale:childrenyes:religiousness ## yearsmarried:childrenyes:religiousness ## gendermale:yearsmarried:education ## gendermale:childrenyes:education ## yearsmarried:childrenyes:education ## gendermale:religiousness:education ## yearsmarried:religiousness:education ## childrenyes:religiousness:education . ## gendermale:yearsmarried:childrenyes:religiousness ## gendermale:yearsmarried:childrenyes:education ## gendermale:yearsmarried:religiousness:education ## gendermale:childrenyes:religiousness:education . ## yearsmarried:childrenyes:religiousness:education ## gendermale:yearsmarried:childrenyes:religiousness:education ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for quasipoisson family taken to be 6.52) ## ## Null deviance: 2925.5 on 600 degrees of freedom ## Residual deviance: 2378.1 on 569 degrees of freedom ## AIC: NA ## ## Number of Fisher Scoring iterations: 9
step(Mod1) #Для квази-пуассоновских моделей эта функция работать # не будет, так как не определен AIC
Mod_nb <- glm.nb(affairs ~ gender * yearsmarried * children *
religiousness + education, data = af)
Anova(Mod_nb)
## Analysis of Deviance Table (Type II tests) ## ## Response: affairs ## LR Chisq Df Pr(>Chisq) ## gender 0.12 1 0.7282 ## yearsmarried 17.87 1 0.000024 *** ## children 1.24 1 0.2647 ## religiousness 15.06 1 0.0001 *** ## education 1.93 1 0.1645 ## gender:yearsmarried 0.22 1 0.6366 ## gender:children 0.02 1 0.8975 ## yearsmarried:children 7.98 1 0.0047 ** ## gender:religiousness 0.01 1 0.9179 ## yearsmarried:religiousness 0.28 1 0.5965 ## children:religiousness 0.02 1 0.8758 ## gender:yearsmarried:children 0.82 1 0.3662 ## gender:yearsmarried:religiousness 0.80 1 0.3715 ## gender:children:religiousness 2.25 1 0.1337 ## yearsmarried:children:religiousness 1.78 1 0.1825 ## gender:yearsmarried:children:religiousness 0.06 1 0.8032 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
step(Mod_nb)
Mod_nb_final <- glm.nb(formula = affairs ~ yearsmarried +
children + religiousness + yearsmarried:children,
data = af, init.theta = 0.1346363532, link = log)
summary(Mod_nb_final)
## ## Call: ## glm.nb(formula = affairs ~ yearsmarried + children + religiousness + ## yearsmarried:children, data = af, init.theta = 0.1346359871, ## link = log) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -1.090 -0.824 -0.700 -0.372 1.910 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -0.4047 0.4309 -0.94 0.3476 ## yearsmarried 0.2978 0.0631 4.72 0.0000024 *** ## childrenyes 1.3847 0.4621 3.00 0.0027 ** ## religiousness -0.4171 0.1064 -3.92 0.0000878 *** ## yearsmarried:childrenyes -0.2232 0.0690 -3.24 0.0012 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for Negative Binomial(0.135) family taken to be 1) ## ## Null deviance: 375.76 on 600 degrees of freedom ## Residual deviance: 335.57 on 596 degrees of freedom ## AIC: 1477 ## ## Number of Fisher Scoring iterations: 1 ## ## ## Theta: 0.1346 ## Std. Err.: 0.0148 ## ## 2 x log-likelihood: -1465.2630
Mod_nb_test <- data.frame(.fitted = predict(Mod_nb_final),
.pears_resid = residuals(Mod_nb, type = "pearson"))
ggplot(Mod_nb_test, aes(x=.fitted, y = .pears_resid)) + geom_point() +
geom_smooth(se = FALSE)
Resid_Mod <- resid(Mod_nb_final, type = "pearson") N <- nrow(af) p <- length(coef(Mod_nb_final)) df <- (N - p) fi <- sum(Resid_Mod^2) /df fi
## [1] 0.723
MyData <- expand.grid(yearsmarried = seq(min(af$yearsmarried),
max(af$yearsmarried)), children = c("yes", "no"),
religiousness = seq(min(af$religiousness), max(af$religiousness)))
MyData$Predicted <- predict(Mod_nb_final, newdata = MyData,
type = "response")
ggplot(MyData, aes(x = yearsmarried, y = Predicted,
color = religiousness)) + geom_line(aes(group = religiousness),
size = 1.5) + facet_grid(~children, labeller = label_both) +
scale_color_gradient(low = "yellow", high = "red") +
geom_point(data = af, aes(x = yearsmarried, y = affairs),
position = position_jitter(width = 1, height = 1))